---

Faculté de génie

Département de génie chimique et de génie biotechnologique

---

Pierre Proulx

ing. Ph.D., professeur

Manipulations Mathematiques Assistées par Ordinateur

Logiciels commerciaux

• Maple

• Mathematica

• Mathcad

Logiciels gratuits, libres, etc...

• Maxima (originalement macsyma, et MAC dans les années 60 au MIT, très solide code avec une syntaxe et un language propre)

• Sage (développé sur python et d'autres sources, maintenant un peu sur le déclin)

• sympy (package python)

• scipy (package python incluant sympy, numpy, etc...)

Pourquoi utiliser sympy et scipy ?

1. Le language python est, selon moi, le language de programmation le plus naturel et intuitif que j'aie vu en 40 ans comme utilisateur, et non comme programmeur.

2. Le language python est largement classé parmi les languages les plus utilisés et les plus enseignés dans les grandes universités américaines depuis plusieurs années.

3. Sympy et scipy fournissent des outils simplifiant de beaucoup les tâches lourdes, répétitives et fastidieuses associées aux manipulations algébriques allant jusqu'à la solution d'intégrales et d'équations différentielles.

---

Que fait sympy? Du calcul symbolique, ce qui veut dire qu'il fait les manipulations avec les variables algébriques tel qu'on le fait nous même. Il ne fait pas les calculs, si on ne lui demande pas.

---

Les quelques pages qui suivent vous donneront des outils pour vous permettre d'alléger les manipulations mathématiques dans tous les cours que vous suivrez pendant votre formation. Il ne tient qu'à vous d'utiliser ces outils qui ne visent pas à remplacer la compréhension mais bien à vous permettre de passer plus de temps à l'application des résultats qu'aux manipulations fastidieuses.

---

---

Algèbre élémentaire

import sympy as sp
from IPython.display import * # pour utiliser display, fonction un peu similaire à print
sp.init_printing(use_latex=True)  # pour l'affichage professionnel de LaTeX.

x=sp.symbols('x')
display(sp.sin(x))
display(sp.sin(x).subs(x,1))
display(sp.sin(x).subs(x,1).evalf(5))

$$\sin{\left (x \right )}$$

$$\sin{\left (1 \right )}$$

$$0.84147$$

Évidemment, sympy fait des tâches comme la recherche de racines de polynômes: $x^3+2x^2+4x+8$

poly=x**3+2*x**2+4*x+8
racines=sp.solve(poly,x)
display(racines)
print(racines) # vous voyez la différence d'affichage?
#
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize']=10,8
# on veut des graphes d'une certaine grandeur
sp.plot((poly,(x,-3,3)));
# le ; enlève  <sympy.plotting.plot.Plot at 0x1155ae198> après le graphe. Testez-le

$$\left [ -2, \quad - 2 i, \quad 2 i\right ]$$

[-2, -2*I, 2*I]

Isoler une variable dans une expression (si c'est possible), donc 1 équation 1 inconnue:

# exemple tiré de Transport Phenomena, chapitre 6
# trouver la vitesse de chute libre d'une petite particule
g,D,rho,rho_P,mu,v,f1,f2=sp.symbols('g,D,rho,rho_P,mu,v,f1,f2')
Re=rho*v*D/mu
f1=4/3*g*D/v**2*(rho_P-rho)/rho
f2=16/Re   # on veut maintenant isoler v dans l'équation f1=f2
#
equation=sp.Eq(f1-f2)
vitesse=sp.solve(equation,v,dict=True)  # solve retourne un 'tuple' avec les solutions
#
# ici on pourrait solutionner pour plusieurs des variables, alors il faut lui demander un dictionnaire
# si on veut avoir le dictionnaire, sinon il ne donnera qu'un tuple
#
display(vitesse[0])    # le premier élément du 'tuple' qui ne contient qu'un élément!
display(vitesse[0][v])    # addresse l'élément du dictionnaire

$$\left \{ v : - \frac{0.0833333333333333 g}{\mu} D^{2} \left(\rho - \rho_{P}\right)\right \}$$

$$- \frac{0.0833333333333333 g}{\mu} D^{2} \left(\rho - \rho_{P}\right)$$

Qu'est-ce que c'est un dictionnaire et à quoi c'est utile?

Imaginez que vous avez une liste de densités d'éléments chimiques:

noms=['eau','benzene','toluene']
rho=[1000,800,900]    # valeurs bidon...
### et ainsi on pourrait utiliser une boucle
for nom in range(0,len(noms)):
    print(noms[nom],rho[nom])

eau 1000
benzene 800
toluene 900

Avec un dictionnaire, on travaillerait directement avec le nom de la substance, pas avec un numéro:

rho={'eau':1000,'benzene':900,'toluene':800}
print(rho)
print(rho['benzene'])

{'eau': 1000, 'benzene': 900, 'toluene': 800}
900

Solution 4 équations a 4 inconnues, pas plus difficile avec sympy, faites-le à la main.

x,y,z,w=sp.symbols('x y z w')
#
systeme=[sp.Eq(x + 5*y +z +2*w,2),
       sp.Eq(-3*x + 6*y -z +15*w,4),
       sp.Eq(-x + 3*y +4*z -10*w,-6),
       sp.Eq(3*x + 16*y -4*z+w,-10)]
print('système')
for i in range(len(systeme)):
    display(systeme[i])    
#
variables=[x, y, z,w]
solution=sp.solve(systeme, variables) # cette ligne solutionne Ax=b
# 
display(solution) # vous voyez ici, par rapport à la solution précédente 
print('solution')
for variable in solution:
    print(variable,'=',solution[variable])

système

$$2 w + x + 5 y + z = 2$$

$$15 w - 3 x + 6 y - z = 4$$

$$- 10 w - x + 3 y + 4 z = -6$$

$$w + 3 x + 16 y - 4 z = -10$$

$$\left \{ w : \frac{1315}{1459}, \quad x : \frac{2197}{1459}, \quad y : - \frac{837}{1459}, \quad z : \frac{2276}{1459}\right \}$$

solution
x = 2197/1459
y = -837/1459
z = 2276/1459
w = 1315/1459

Approximations en séries, tracer des graphes de fonctions symboliques

# Série de Taylor
display('Série de Taylor (8eme ordre) de la fonction erreur',
        sp.erf(x).series(x,0,8))
# On utilise le removeO() qui enlève le terme d'erreur
display(sp.erf(x).series(x,0,8).removeO())

'Série de Taylor (8eme ordre) de la fonction erreur'

$$\frac{2 x}{\sqrt{\pi}} - \frac{2 x^{3}}{3 \sqrt{\pi}} + \frac{x^{5}}{5 \sqrt{\pi}} - \frac{x^{7}}{21 \sqrt{\pi}} + \mathcal{O}\left(x^{8}\right)$$

$$- \frac{x^{7}}{21 \sqrt{\pi}} + \frac{x^{5}}{5 \sqrt{\pi}} - \frac{2 x^{3}}{3 \sqrt{\pi}} + \frac{2 x}{\sqrt{\pi}}$$

taylor4=sp.erf(x).series(x,0,4).removeO()
taylor6=sp.erf(x).series(x,0,6).removeO()
taylor16=sp.erf(x).series(x,0,16).removeO()
#
p=sp.plot((sp.erf(x),(x,0,2)),
        (taylor4,(x,0,2)),
        (taylor6,(x,0,2)),
        (taylor16,(x,0,2)),ylabel='erf(x)',legend=True,show=False)
p[0].label='erf'
p[0].line_color='red'
p[1].label='Taylor ordre 4'
p[1].line_color='blue'
p[2].label='Taylor ordre 6'
p[2].line_color='black'
p[3].label='Taylor ordre 16'
p[3].line_color='green'
p.show()

---

Calcul différentiel et intégral

f=sp.cos(x)
display(f)

$$\cos{\left (x \right )}$$

display(sp.Derivative(f,x))    # affiche la dérivation par rapport à x

$$\frac{d}{d x} \cos{\left (x \right )}$$

display(sp.diff(f,x))          # effectue la dérivée par rapport à x

$$- \sin{\left (x \right )}$$

display(sp.diff(f,x).subs(x,1))          # effectue la dérivée par rapport à x en un point

$$- \sin{\left (1 \right )}$$

display(sp.diff(f,x).subs(x,1).evalf(4)) # effectue la dérivée par rapport à x en un point et évalue numériquement à 4 chiffres

$$-0.8415$$

display(sp.Integral(f,x))      # affiche l'intégrale par rapport à x

$$\int \cos{\left (x \right )}\, dx$$

display(sp.Integral(f,(x,sp.pi,sp.oo)))      # affiche l'intégrale par rapport à x avec les bornes pi et infini

$$\int_{\pi}^{\infty} \cos{\left (x \right )}\, dx$$

display(sp.integrate(f,x))     # effectue l'intégration indéfinie par rapport à x

$$\sin{\left (x \right )}$$

display(sp.integrate(f,(x,1,2)))  # effectue l'intégration définie par rapport à x entre les bornes

$$- \sin{\left (1 \right )} + \sin{\left (2 \right )}$$

display(sp.integrate(f,(x,1,2)).evalf(4))  # effectue l'intégration définie et évalue numériquement

$$0.06783$$

Le volume d'une sphère est obtenu à partir de la géométrie analytique:

Le volume d'un élément différentiel est:

$dV = r^2 sin(\theta) d\theta d\phi$

que l'on integrera: $\int_{0}^{2 \pi}\int_{0}^{\pi}\int_{0}^{R} r^{2} \sin{\left (\theta \right )}\, dr\, d\theta\, d\phi$

# intégrale triple, un bon exemple, le calcul du volume d'une sphère. 
#
r,theta,phi,R=sp.symbols('r,theta,phi,R')
volume_différentiel=r**2*sp.sin(theta)
display(sp.Integral(volume_différentiel,
                    (r,0,R),
                    (theta,0,sp.pi),
                    (phi,0,2*sp.pi)))
display('Volume ',sp.integrate(volume_différentiel,
                     (r,0,R),
                     (theta,0,sp.pi),
                     (phi,0,2*sp.pi)))

$$\int_{0}^{2 \pi}\int_{0}^{\pi}\int_{0}^{R} r^{2} \sin{\left (\theta \right )}\, dr\, d\theta\, d\phi$$

'Volume '

$$\frac{4 \pi}{3} R^{3}$$

---

Équations différentielles.

Sympy permet de solutionner facilement les équations les plus souvent rencontrées en génie chimique. Pas toutes, mais une grande partie.

---

Équation de Newton, ballistique

### Exemple classique, équation de Newton F=Ma
t,g,v_0,x_0=sp.symbols('t g v_0 x_0')
x=sp.Function('x')(t)
eq=x.diff(t,t)-g         # équation ballistique, trajectoire d'un projectile
display(eq)
x=sp.dsolve(eq,x).rhs           # coté droit de l'équation qui contient l'expression désirée
display(x)
# Pour trouver les deux constantes, on utilise les conditions initiales 
Ci1=x.subs(t,0)-x_0              # en t=0, le projectile est à la position x_0
Ci2=sp.diff(x,t).subs(t,0)-v_0    # en t=0, le projectile est à la vitesse v_0
#
constantes=sp.solve([Ci1,Ci2],sp.symbols('C1,C2'))
#
display(x.subs(constantes))

$$- g + \frac{d^{2}}{d t^{2}} x{\left (t \right )}$$

$$C_{1} + C_{2} t + \frac{g t^{2}}{2}$$

$$\frac{g t^{2}}{2} + t v_{0} + x_{0}$$

---

Faculté de génie

Département de génie chimique et de génie biotechnologique

---

Pierre Proulx

ing. Ph.D., professeur

---

Réacteur chimique isotherme, deux réactions consécutives

$\frac {dC_A}{dt}=-k_1C_A$,$\frac {dC_B} {dt}=k_1C_A-k_2C_B$ et $\frac {dC_C}{dt}=k_2C_B$

t,C_A0,C_B0,C_C0,k_1,k_2=sp.symbols('t,C_A0,C_B0,C_C0,k_1,k_2')
C_A=sp.Function('C_A')(t)
C_B=sp.Function('C_B')(t)
C_C=sp.Function('C_C')(t)
e1=sp.Eq(sp.diff(C_A,t)+k_1*C_A,0)           
e2=sp.Eq(sp.diff(C_B,t)+k_2*C_B-k_1*C_A,0)
e3=sp.Eq(sp.diff(C_C,t)-k_2*C_B,0)
#
display(e1,e2,e3)
#

$$k_{1} \operatorname{C_{A}}{\left (t \right )} + \frac{d}{d t} \operatorname{C_{A}}{\left (t \right )} = 0$$

$$- k_{1} \operatorname{C_{A}}{\left (t \right )} + k_{2} \operatorname{C_{B}}{\left (t \right )} + \frac{d}{d t} \operatorname{C_{B}}{\left (t \right )} = 0$$

$$- k_{2} \operatorname{C_{B}}{\left (t \right )} + \frac{d}{d t} \operatorname{C_{C}}{\left (t \right )} = 0$$

C_A,C_B,C_C=sp.dsolve([e1,e2,e3])      
#
display(C_A,C_B,C_C)

$$\operatorname{C_{A}}{\left (t \right )} = C_{1} e^{- k_{1} t}$$

$$\operatorname{C_{B}}{\left (t \right )} = \frac{C_{1} k_{1} e^{- k_{1} t}}{- k_{1} + k_{2}} + C_{2} e^{- k_{2} t}$$

$$\operatorname{C_{C}}{\left (t \right )} = - \frac{C_{1} k_{2} e^{- k_{1} t}}{- k_{1} + k_{2}} - C_{2} e^{- k_{2} t} + C_{3}$$

#
Ci1=C_A.subs(t,0).rhs-C_A0    # attention, seulement la droite de l'équation nous interesse
Ci2=C_B.subs(t,0).rhs-C_B0
Ci3=C_C.subs(t,0).rhs-C_C0
constantes=sp.solve([Ci1,Ci2,Ci3],sp.symbols('C1,C2,C3'))
display(constantes)

$$\left \{ C_{1} : C_{A0}, \quad C_{2} : \frac{1}{k_{1} - k_{2}} \left(C_{A0} k_{1} + C_{B0} k_{1} - C_{B0} k_{2}\right), \quad C_{3} : C_{A0} + C_{B0} + C_{C0}\right \}$$

C_A=C_A.subs(constantes)
C_B=C_B.subs(constantes)
C_C=C_C.subs(constantes)
display(C_A,C_B,C_C)     # la solution est obtenue

$$\operatorname{C_{A}}{\left (t \right )} = C_{A0} e^{- k_{1} t}$$

$$\operatorname{C_{B}}{\left (t \right )} = \frac{C_{A0} k_{1} e^{- k_{1} t}}{- k_{1} + k_{2}} + \frac{e^{- k_{2} t}}{k_{1} - k_{2}} \left(C_{A0} k_{1} + C_{B0} k_{1} - C_{B0} k_{2}\right)$$

$$\operatorname{C_{C}}{\left (t \right )} = - \frac{C_{A0} k_{2} e^{- k_{1} t}}{- k_{1} + k_{2}} + C_{A0} + C_{B0} + C_{C0} - \frac{e^{- k_{2} t}}{k_{1} - k_{2}} \left(C_{A0} k_{1} + C_{B0} k_{1} - C_{B0} k_{2}\right)$$

# essai avec des valeurs 
dico={'C_A0':1,'C_B0':0.5,'C_C0':0,'k_1':.1,'k_2':0.08}
p=sp.plot((C_A.rhs.subs(dico).evalf(),(t,0,30)),
       (C_B.rhs.subs(dico).evalf(),(t,0,30)),
       (C_C.rhs.subs(dico).evalf(),(t,0,30)),ylabel='C',legend=True,show=False);
p[0].label='C_A'
p[0].line_color='red'
p[1].label='C_B'
p[1].line_color='blue'
p[2].label='C_C'
p[2].line_color='black'
p.show()